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Abstract — We study the problem of optimal estimation 
using quantized innovations, with application to distributed 
estimation over sensor networks. We show that the state 
probability density conditioned on the quantized innovations 
can be expressed as the sum of a Gaussian random vector 
and a certain truncated Gaussian vector. This structure bears 
close resemblance to the full information Kalman filter and so 
allows us to effectively combine the Kalman structure with a 
particle filter to recursively compute the state estimate. We call 
the resuting filter the Kalman like particle filter (KLPF) and 
observe that it delivers close to optimal performance using far 
fewer particles than that of a particle filter directly applied to 
the original problem. We also note that the conditional state 
density follows a, so called, generalized closed skew-normal 
(GCSN) distribution. 

Index Terms — Distributed state estimation. Sign of Innova- 
tion, Closed Skew Normal Distribution, Particle Filter, Wireless 
sensor network, Kalman Filter. 

I. Introduction 

Recent advances in very large-scale integration and mi- 
croelectromechanical system technology have led to the 
availability of cheap, low quality and low power consumption 
sensors in the market. This has generated a great deal of 
interest in wireless sensor networks (WSNs) due to their 
potential applications in several diverse fields [1]. Sensor 
network constraints such as limited bandwidth and power 
have inspired a considerable amount of research in de- 
veloping energy efficient algorithms for network coverage 
and decentralized detection and estimation using quantized 
sensor observations [2]-[4]. 

The problem of estimation with quantized measurements 
is almost as old as the Kalman filter itself. An early survey 
on the subject can be found in [5]. However, most of the 
earlier techniques centered around using numerical integra- 
tion methods to approximate the optimal state estimate. The 
advent of particle filtering [6]-[8] created a whole new set of 
tools to handle non-linear estimation problems. For example, 
[4] proposes a particle filtering solution for optimal filtering 
using quantized sensor measurements. But, quantizing sensor 
measurements can lead to large quantization noises when the 
observed values are large which then leads to poor estimation 
accuracy. In [9], this limitation is overcome by developing an 
elegant distributed estimation approach based on quantizing 
the innovation to one bit (the so called sign of innovation 
or SOI). In [10], this is generalized to handle multiple 
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quantization levels. In both cases, it is assumed that the 
conditional state density is approximately Gaussian leading 
to a linear filter and a very simple characterization of its 
error performance. Under the Gaussian assumption, the error 
covariance matrix associated with the state estimation error 
satisfies a modified Riccati recursion of the type that appears 
in [11]. The only difference between this modified Riccati 
and the traditional one is a scaling factor A multiplying 
the nonlinear term of the recursion. For the SOI Kalman 
filter (SOI-KF), X is ^ while [10] presents a formula for 
A in the case of multiple quantization levels. Henceforth, 
these filters will be referred to as SOI-KF and MLQ-KF, 
and their associated Riccati recursions as SOI-Riccati and 
MLQ-Riccati respectively. 

For linear time invariant dynamical systems, if the Gaus- 
sian assumption were realistic, convergence of the modified 
Riccati must mean the convergence of the corresponding 
linear filters. Using results presented in [11] one can find 
linear time invariant systems for which the MLQ-Riccati 
and SOI-Riccati converge. [12] provides examples for which 
the actual error performance of SOI-KF and MLQ-KF do 
not converge to their respective Riccatis, which warrants a 
closer examination of the assumption of Gaussianity. This 
is one of the questions addressed in this paper. We derive 
a novel stochastic characterization of the conditional state 
density (see theorem [3. lb . Using this, it is straighforward to 
see that the conditional state density is not Gaussian, as one 
would expect, given the non-linear nature of quantization. In 
fact, it belongs to a class of distributions which we refer to 
as Generalized Closed Skew Normal (GCSN) distributions. 
A careful literature review reveals that a related observation 
has been made in [13]. In particular, with some effort, [13] 
can be used to derive theorem 13.11 Specializing this result 
to state space models, we develop a novel particle filtering 
approach to optimally estimate the state using quantized 
measurements/innovations. In the rest of the paper, we use 
the words 'measurements' and 'innovations' interchangealy 
since the analysis will prove that the general structure of the 
filter does not depend on whether sensor measurements or 
innovations are quantized. The proposed filter requires far 
fewer particles than that of a particle filter applied directly 
to the original problem [12], as will be shown through 
various simulations. We also develop a precise formulation 
of the conditional state density and observe that it follows 
what we call a generalized closed skew-normal distribution, 
which is very similar to those studied in [14]-[19]. Some 
useful properties of this distribution are also provided in the 
Appendix. The next section introduces the problem setup. 
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Fig. 1. WSN with a fusion center; The sensors act as data gathering 
devices. 



II. Problem Statement and Preliminaries 

The broader problem that one would like to solve can 
be cast as causal estimation of a random process {x{n)} 
using a quantized version of a measurement process {y{n)}, 
where x{n) and y{n) can be vectors. Without any further 
structure, this is a difficult problem to analyze. When {x{n)} 
and {y{n)} are jointly Gaussian, we will provide a novel 
stochastic characterization of the probaility density of x{n) 
causally conditioned on the quantized measurement process. 
We use it to propose a novel filtering technique for the above 
problem which reduces to an elegant particle filter when 
{x{n)} has state space structure. 

A. Notation 

The following notation will be used in the rest of the paper. 

1) If {w(«)}55^_oo is a discrete time random process, 
u{i : j) denotes the collection of random vari- 
ables (u(i) , . . . , u{j)) . Un denotes the random vector 
[u(0), . . . , u{n)]^ and u„ denotes a realization of Un- 

2) For a random vector y, \\Yf = E{Y -EY){Y -EY'f 

3) For random vectors X, Y, {X, Y) = E{X - EX){Y - 

eyY 

4) For random variables {Xi, . . . , Xn), C(Xi, . . . , X„) 
denotes their linear span. 

5) For a vector x, x^*'' denotes its i-th component. In the 
context of particle filtering, a:* would denote the z-th 
particle. 

6) Nkii-J., S) denotes a k-dim normal random variable with 
mean /i and covariance E. Nk{a,b, fj.,Y.) denotes a k- 
dim normal truncated to lie in (a, b), where a, b are 
fc-dim vectors and the truncation is component-wise. 

7) ^{x) = P{X < x), where X - iV(0, 1), 
$(a,6,/i,cr) = P{X e (a, 6)) when X - N{^,a^) 
and in general, $(J*, /i, cr) = P{X e I*). 

8) The notion of optimality to be used throughout the paper 
is mean squared error optimality. 

B. Problem Formulation 

Suppose {x{n)} and {y{n)} have the following state space 
structure. 



where x[n) G is the state, y[n) G M is the observation, 
and w{n) G M'' and v{n) G M are uncorrelated Gaussian 
white noises with zero means and covariances W{n) and 
R{n) = (j1{'n), respectively. The initial state, a;(0), of the 
system, is also a zero mean Gaussian with covariance P(0) 
and is uncorrelated with both w{ri) and v{n). As in [9], [10], 
we consider the sensor network configuration in which the 
fusion center has sufficient power to broadcast its predicted 
output and the corresponding error covariance to its sensors. 
Sensors are assumed to have limited power and hence their 
transmission of information should be limited. Here, we 
assume that the energy required for receiving messages is 
much less than that for transmitting. 

Fig [T] outlines the overall filtering paradigm. Once a 
scheduling algorithm is in place, at each time instant, a 
sensor S{n) makes a measurement y{n) and computes the 
innovation e(n) = y{n) — y{n\n — 1), where y{n\n — 1) = 
hx{n\n — 1) together with the variance of the innovation 
||e(n)|p are received by the sensor from the fusion center, 
with x{n\n—l) being the one step predictor of the state. [9], 
[10] propose methods to quantize e{n) and use the quantized 
innovations to update the state estimate. These filters take the 
following shape 

, , , , , , , , Pln){H{n)Y 
xinin) ^ xinin - 1) + L (g(n)) Hin)P{n)iH{n))T + 

x{n + 1/n) = A{n)x{n / n) 

P{n){H{n)Y H{n)P{n) 



P{n/n) = P{n) - X- 



(2a) 



x{n + 1) = A{n)x[n) + 'w{n) 
y{n) = H(n)x{n) + v(n) 



(la) 
(lb) 



Hin)Pin)iHin))T + cT^ 

P{n + 1) = P{n + 1/n) = A{n)P{n/n){A{n)f + W{n) 

(2b) 

where q{n) denotes the quantized innovation while L {q{n)) 
and the value of A depend on the quantization scheme used. 
Eqs ( l2al i and (IZbl i constitute the modified Riccati recursion 
with parameter A. For SIO-KF, A = | and for MLQ-KF, 
[10] provides a formula for A and L{q{n)). The above 
filter is derived based on the assumption that the conditional 
distribution / (x(n)/q„_i) is Gaussian, which we will prove 
is generally false. [12] provides examples where the error 
performance of the filters in [9], [10] do not track the modi- 
fied riccati recursions that they were predicted to. Instead, the 
optimal filter, which was approximated by a particle filter, 
has been observed to have an error covariance matrix that 
obeys the modified Riccatis. In order to approximate the 
optimal filter, [12] employs a very basic particle filtering 
algorithm which is outlined below for easy reference. 

Alg 1. Particle Filter 

1) Set n = 0. For i = 1, ■ • • ,N, initialize the particles, 
x' (0| - 1) ~ ./(a;(0)) and set x{0\ - 1) = 

2) At time n, set q{n) = Q {y{n) — H {n)x{n\n — 1)), 
where Q{.) is a quantization rule. 

3) Suppose the quantized value q{n) implies that 
y{n) — H {n)x{n\n — 1) G liii), then 

v{n) + H{n){x{n) — x{n/n — 1)) G 2{n). The 



importance weights are now calculated as follows 

w\n) = <^{I{n),H{n){x{n) - x{n/n - 1)), cT„(n)) 
4) Measurement update is given by 



N 

x{n/n) = — 1) 

where W''{n) are the normalized weights, i.e., 

Zjlli™ (") 

5) Resample iV particles with replacement accoding to, 

Prob {x^{n/n) ~ [n/n — 1)) = [n) 

6) For i = 1, • • • , A^, predict new particles according to, 

x\n + 1/n) ~ / {x(n + l)|a;*(n|n)) , i.e., 
x'{n + 1/n) = A{n)x\n/n) + iVfc(0, W{n)) 

where W{n) is the covariance of the process noise at 
time n. 

7) Set x{n + \ \n) ~ A{n)x{n\n) . Also, set n = n + 1 
and iterate from step 2. 

The particles in Alg 1 describe the conditional state 
density / (a;(n)/q„) and simulations suggest that one needs 
upwards of a few hundred particles to get satisfactory error 
performance for most systems. In the following sections, 
we disprove the premise behind MLQ-KF and SOI-KF and 
develop a novel particle filtering technique (KLPF) which 
converges to the optimal filter asymptotically. Simulations 
show that the KLPF needs far fewer particles than the particle 
filter of Alg 1 . The difference partly lies in using particles to 
describe a probability density funtion with much less support 
than the conditional state density. 

III. Full Information Vs Quantized Innovations 

Suppose {x{n)} and {y{n)} are jointly Gaussian, then it is 
well known that the probability density of x{n) conditioned 
on y.n is a Gaussian with the following parameters 

x{n)/ - y„) ^ Nk{Q, Ra) + Rxvin){Ry{n))-\^, 

(3) 

and Ra ^ Rx{n) - Rxyin){Ry{n))~'^ Ryx{n) 

When {x{n)} has an underlying state space structure and 
{y{n)} is a linear measurement of {x{n)} corrupted by 
additive white Gaussian noise, as defined in Eq ([T]i, it is 
well known that the following Riccati recursion propogates 
the error covariance P{n) = Ra 

P{n) = P{n/n - 1) = A{n)P{n -l/n- l)A(n)^ + Q(n) 
P{n)H{n){H{n))^P{n) 



P{n/n) = P{n) 
P(0) = ||x(0)f 



H{n)P{n){H{n)y 



(4) 



(5(.) be a general quantization rule which generates the quan- 
tized measurement process {q{n)}, where q{n) = Q(yn)- 
Note that we allow the quantization rule to depend on all 
the past measurements y{i), i < n. This includes, as a 
special case, the method of quantizing the innovations first 
proposed in [9]. We will show that the probability density of 
x{n) conditioned on the quantized measurements Q„ admits 
a characterization very similar to Eq (|3]l. We state the result 
in the following theorem. 

Theorem 3.1: The probability density of x{n) conditioned 
on the quantized measurement Q„ can be expressed as 

x{n)/{Qn^qn) -Nk{0,RA) 

+ Rxy{n){Ry[n))-^ (3^n/ (Q„ = q„)) 

(5) 

Proof: The proof is fairly straightforward. The theorem 
will be proved by showing that the moment generating 
fuction of x{n) / [Qn — cin) can be seen as the product 
of two moment generating functions corresponding to the 
two random variables in Eq (|5]l. For brevity, we will write 
x{n)l (Q„ = q„) as x(n)/q„. 



/(x(n)/q„)=y /(.T(n),y„/q„)dy„ 
Noting that /(a:(n)/y„, q„) = /(a;(n)/y„), we can write 
i;e*"-("Vqn = / e*"-(")/(x(n)/y„)/(y„/q„)d.T(n)dy„ 



'/(yn/qn)rfy n 



We would like to address the problem of optimal estimation 
using a quantized version of the observation process. Let 



^ mfg of i?.„(n)(i?„(n))-iyn/q„ 

=^ M,(^„y^^{t) = Mz{t)Myj^^{{Ry{n))-^Ryx{n)t) 

(6) 

where Z ^ Nk{Q, Ra)- In getting (*), we used the fact that 
x{n)/yn ~ Nk{Rxy{n){Ry{n))~^yn, Ra)- For any random 
variable Y, it is easy to see that MyiA^t) = MAvit)- The 
result is now obvious from Eq (|6]l. ■ 

Comparing Eqs ^ and Q, the only difference is the mea- 
surement vector y„ being replaced by the random variable 
yq{n) = yn/ln- It IS casy to see that 3^q(n) is a multivariate 
gaussian random variable truncated to lie in the region 
defined by q„. It is worth noting that the covariance of 
x{n)/qn, ||a;(r7,)/q„p, is given by 

Ra + Rxy{n){Ry{n))-^\\yr,/qnf{Ry{n))-^Ryx{n) 

As the quantization scheme becomes finer (in an appropriate 
manner), 3^q(n) clearly converges to Xi and a;(n)/q„ ap- 
proaches a Gaussian as is well known. Using theorem l3.ll 
it is easy to see that 2;(n)/q„ is not Gaussian in general, 
contrary to the assumption made in [9], [10]. It belongs 
to a class of distributions, which we call the Generalized 
Skew Normal Distributions (GCSN), the details of which 
are provided in the Appendix. 



We will use theorem 13.11 to propose a novel particle 
filtering scheme. We begin by noting that 



Ex{n)/cin = R^y{n){Ry{n))~^Ey^{n) 



(7) 



The filter is implemented for the quantization scheme pro- 
posed in [9], [10]. At time n, the sensor which is sched- 
uled to make the measurement receives a prediction of its 
measurement y{n) = y{n/n — 1) and the error covariance 
\\y{n) — y{n/n — The sensor, then quantizes pl^i 

where e(n) = y{n) — y{n/n — 1), using the following 
quantizer and transmits the result to the fusion center. 



Q{x) 



if a; > rx-i 

if rj_i < X < ri, i >2 

if X < ri 



where (ri, . . . , tk-i) are the quantization levels. Using the 
convention, ro = — oo and rx = cw, we can re-write 
the quantization levels as (rp, . . . , ). The output of the 
quantizer at time n will be denoted by q{n) and the lower and 
upper limits of the interval implied by q{n) will be denoted 
by ri (n) and r2 (n) respectively. We assume that the fusion 
center has access to the exact value of q{n) at every instant. 
Note that y{n/n — 1) is an estimate of the measurement, not 
necessarily the optimal. Hence y{n) — y{n/n — 1) is not an 
innovation in the true sense of the word, unless the estimator 
employed at the fusion center is optimal. Nevetheless, we 
will refer to q{n) as the quantized innovation. With this 
setup, we will develop a filtering technique which needs 
far fewer particles to achieve optimal performance than the 
simple particle filter applied to the original problem. 

IV. The Filter 

We shall begin with the observation that 3^q(n) is a 
multivariate normal distribution truncated as follows 

= yn/{yij) e (ri(j) + r2(j) + V j < 

Let si{n) ^ ri{j) + y{j) and S2{n) ^ r2{j) + y{j). 
Then, clearly, 3^q(n) is a multivariate normal with its j-\h 
component truncated to lie in the interval (si(j), S2(j)) V 
j < n, i.e., yq{n) - 7V„+i(si,„, S2,„, 0, i?2/(n)). From Eq 
O, it is clear that the optimal state estimate can be computed 
by first computing the mean of the truncated normal. Before 
proposing the filter, we need a couple of results for 

(a) relating the distributions of yqin) and 3^q(ri — 1) and 

(b) generating scalar truncated normal random variables. 
Lemma 4.1: Let (Zi, Z2, • • • , Z„) ~ /„ (z„), where 

fn (Zri) = Nn (z„ ; Si_„ , S2,„ , 0, i?^ (n) ) , 

then the marginal of the first n — 1 components is given by 

{Zi, • • • , Zn-l) ^ fn (Zri-l) 
=/„-l(z„-i) 



$ {si{n), S2{n),Ez{n)/z„-i,y/RAn^ 
Also the conditional distribution of Z„/Zi:„_i is given by 

/„ (z(n)/z„_i) =iV(^z(ri); Si(n), S2{n), Ez{n) /zn-i, Rau 

Ez{n)/zn-i = {z{n), Zn-i){Rz{n-l)y^Zn-i 

RAn = Mn)f - (^(n),Z„_i)(i?,(7i-l))-i(Z„_i,z(n)) 

Note that /„ (z(7i)/z„_i) is a one dimensional truncated 
normal. 

Proof: The proof is straightforward and is not presented 
here due to space limitations. ■ 

The following Lemma describes a standard technique to 
generate scalar truncated normal random variables. 

Lemma 4.2: Let U ^ U{0, 1), then 

Y ^ (($(6) - $(a)) U + $(a)) 

is distributed as N{a, b, 0, 1) 

Proof: If F{y) denotes the cumulative distribution 
function of Y, then it is well known that F{Y) ~ U{0, 1). 
We also have 



if ?; < a 



1 



if y>b 



So, if [/ ~ C/(0,1), then F^\U) - iV(a,5,0, 1) 
and hence the technique. Noting that N{a,h, ^,a'^) = 



N{ 



a — fi b—fi 



, 0, 1), we can generate scalar truncated normal 



random variables with arbitrary mean and variance. ■ 
We will now propose a particle filtering technique to 
recursively compute the state estimate. 



Alg 2. Truncated Normal Particle Filter 

1) At n ~ 0, generate 

{y'mti - iV(si(0),S2(0),0,i?,y(0)) using the 
technique outlined in Lemma W2\ 

2) At time n, for each particle {yJi_i}, compute the 
weight as 

w'{n) = $ (si(n), S2(n), Ey{n)/yl_-^^, y/R^^^ (8) 

3) Generate 

{y\n)}l, ^ N {si{n),S2{n),Ey{n)/rn_„RAn) 
and define 



4) Measurement update 



x{n\n) = R.,y{n){Ry{n))-'^^^^^^ (9) 

5) Resample the N particles {yj^l^i with replacement 
accoding to Prob (yjj = where the normalized 

weights are given by 



CC Nn-l (Z 0,i?,(n-l)) 



6) Set x{n + l\n) = A{n)x{n\Ti). Also, set n = n ■ 
and iterate from step 2. 



In the filter proposed above, it is important to note that the 
dimension of the particle y*j is n and hence increases with 
time. In the absence of any structure on the random processes 
{x{n)} and {2/(71)}, the computations involved in running 
the above filter will quickly become infeasible. Alg 2 is of 
purely theoretical interest but the technique presented above 
is greatly simplified when x{n) has state space structure. In 
the following subsection, we will show how to overcome 
the problem of increasing particle dimension when {x{n)} 
is a Gauss Markov process and {y{n)} is a linear mea- 
surement of {x{n)} corrupted by additive white Gaussian 
noise. The filter takes an elegant formulation in which the 
particle dimension remains fixed and is equal to the state 
dimension. This approach requires far fewer particles than 
Alg 1 and the filtering technique is quite general, in that, 
it can handle arbitrary number of quantization levels and 
arbitrary quantization intervals. We will also observe that 
the filter requires fewer and fewer particles as the number of 
quantization levels increases. 



A. Exploiting State Space Structure - The Kalman Like 
Particle Filter 

Suppose {x{ny\ and {y{ny} have the state space structure 
defined in Eq ([TJ. Then we know that Rxy{n){Ry{n))~^yl^ 
is the optimal estimate of x{n) given the measurement 
vector and hence can be computed by running a Kalman 
filter using y^, the specific realization of the measurement 
vector y„. Similarly Ey{n)/yl^_i is the optimal estimate 
of y(n) given 3^„_i = y,\_i and Rau is the resulting 
error covariance. All these quantities emerge naturally in the 
following Kalman filtering steps. 

..,0/0, = .^^^^^'T'"';, 



Now consider Eq (|9]l, we can write it alternately as 



H{0)PmH{0)V + aUOr 
P{k + 1) = P{k + 1/fc) = AP{k/k)A^ + W{k) 



(10a) 
(10b) 



P(fc + l/fc + 1) = P{k) - 



P{k)H{k){H{k)Y P{k) 
H{k)P{k){H{k)Y + al{k) 



x\k + 1) ^ x\k + 1/fc) = Ax^{k/k) 
x^k + l/k+l) = x\k+l)+ 
P{k){H{k)r 



(10c) 
(lOd) 
(lOe) 



H{k)P{k){H{k)Y + <jl{k) 



{y\k + l)-H{k + l)x\k+l)) 
(lOf) 



Now, R^y{n){Ry{n)) Vji, Rau and Ey{n)/yi^_-^ can be 
calculated as follows 

R,y{n){Ry{n))-'yl-^ = x\n/n) (11a) 
Rau = H{n)P{n){Hin)f + alin) (lib) 
Eyin)/rn_,^H{n)x\n) (11c) 



x{n\n) = R.^.^(n){Ry{n))-'^^^^^^ 



N 



N 



= Y,w'{n)R,y{n){Ry{n)r^rn = J] W^Hx^n/n) 

i=l 1=1 

(12) 

From Eqs ( fTTT i and ( fT2b . it is easy to see that all the 
information in y,'j is captured in x^injn). Hence, we only 
need to work with {a;*(n/n)}^j^ at any time n. We can now 
desribe the new filter as follows. 

Alg 3. Kalman Like Particle Filter (KLPF) 

1) At n = 0, generate 

{2/*(0)}iIi - iV(si(0),S2(0),0,i?,(0)) using the 
technique outlined in Lemma |4!2] Compute {a;*(0)} 
using Eq (fTOal i and x*(l/0) = A(0)2;'(0/0) 

2) Generate 

{y'Hiili ~ N (si(n), S2(7i), H{n)x\n), Rau) 



Use Eq ( II Ibl l to compute i?ATi- Then form {x'^{n/n)} 
using Eq (fTol l. 

3) At time n, for each particle {a;*(ri)}, compute the 
weight {w*(ri)} as 

w'{n) =$(^si(ri),S2(ri),i?(n)a;*(n), i/^An) 

4) {Measurement update) 
Normalize the weights to get iZ;'(n) 



and 



compute the measurement updated state estimate using 
Eq ( fT2l ). i.e., x{n/n) = '^^^iw'' [n)x^ [n/n) 

5) Resample the N particles {x*(n/n)}^j^ with 
replacement accoding to Prob i^x'^in/n)) = vifin) and 
compute a-*(n + 1) = A(n + \)x'^{n/n). 

6) Set x{n + l|n) = A{n)x{n\n). Also, set n = n + 1 
and iterate from step 2. 

Remark 1: It is easy to see that the above filter can be 
extended to the case when the quantizer employed at the 
sensor is of the more general form Q{x) ~ qj if x £ Xj for 
1 < < K- 

From the above implementation, in terms of complexity, the 
KLPF is clearly equivalent to running N parallel Kalman 
filters. Hence, the complexity of the KLPF scales linearly 
in N , the number of particles. Also, it converges to the 
optimal filter as iV ^ 00. But for most systems, simulations 
suggest that the KLPF delivers close to optimal performance 
for < 50. The particles in the KLPF describe the random 
variable Rxy{n-){Ryin))'^yq{n). The support of its distri- 
bution clearly decreases with increasing quantization levels. 
As a result KLPF needs fewer particles as the quantization 
becomes finer, a property that Alg 1 does not share. This 
will be demonstrated through examples in Section [V] 

The scenario considered thus far involves one measure- 
ment per time instant. But Alg 3 can be easily extended to 
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Summary of numerical results 
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Example 2 


SOI 


2 bit 


SOI 


2 bit 


SOI-KF 


no 




yes 




MLQ-KF 




no 




yes 


Alg 1 


2500 


10000 


500 


750 


KLPF 


500 


90 


25 


3 



-Monte Carlo of SOI - KF 
-SOI-Riccati 
.Particle Filter, N„p = 2500 




(a) SOI-KF, Alg 1 and KLPF 



handle multiple measurements from different sensors at a 
given time instant. Here it is assumed that the measurement 
noise processes are uncorrected across different sensors. 
The proof is simple and is not detailed here due to space 
constraints. 

V. Simulations 

In Alg 1, the particles describe the full probability density 
of the state conditioned on quantized measurements. While 
in Alg 3, part of the information about the conditional state 
density is captured neatly by the Kalman filter. So, the 
particles describe a truncated Gaussian which has much less 
support than the full conditional distribution. We give two 
examples in this section to demonstrate the effectiveness of 
KLPF. The following table summarizes the highlights from 
the two examples 

In Table U a 'yes' indicates that the filter works and is 
close to optimal, a 'no' indicates that its estimation error 
diverges and a '-' means that the quantization method does 
not apply to the filter 'SOI' stands for 'sign of innovation' 
and '2-bit' stands for a quantization rule with quantization 
intervals given by (-oo, -1.2437), (-1.2437,-0.3823), 
(-0.3823,0.3823), (0.3823,1.2437) and (1.2437, oo). If 
the innovation falls in the interval (—0.3823, 0.3823), no 
measurement update is done, so that 2 bits will suffice to 
represent the output of the above quantizer. The numbers in 
front of Alg 1 and KLPF denote the number of particles 
required to approximate the optimal filter Clearly, KLPF 
requires far fewer particles than Alg 1. Also evident from 
Table U is the fact that KLPF needs dramatically fewer 
particles as the quantization becomes finer. 

A. Example 1 

We re-use this example from [12]. Consider a linear 
time invariant system of the form ([T]) with the following 

W = 2Z3, 
denotes an 



0.95 1 
0.9 10 
0.95 



102 



parameters: A 

R ^ al ^ 2.5 and P(0) = O.GII3, where 1,^ 
m X m identity matrix. Note that A is a stable matrix. As 
can be seen from the plots, SOI-KF and MLQ-KF diverge 
but KLPF delivers optimal performance with much fewer 
particles than Alg 1. With the addition of just 1 bit, the 
required number of particles drops from 500 to 90. But with 
Alg 1, the number of particles goes up from 2500 to 10000. 

Remark 2: Extensive simulations suggest that the mod- 
ified Riccati of Eq (O describes the error performance of 



— MLQ-Riccati for 2 bits 

— Monte Carlo Simulation of MLQ-KF 

.Particle Filter, N^^ = 10000 



(b) MLQ-KF, Alg 1 and KLPF 

Fig. 2. In (a), SOI-KF clearly diverges while Alg 1 and KLPF 
converge to the optimal filter From (b), MLQ-KF with 4 levels of 
quantization also diverges while KLPF converges to the optimal 
filter with just 90 particles 



the optimal filter whereas MLQ-KF is itself clearly not the 
optimal filter In Examples 1 and 2, by optimality, we meant 
that the mean squared errors of Alg I and KLPF tracked 
the corresponding modified Riccati recursions very closely. 
In fact, simulations showed that the mean squared errors of 
Alg 1 and KLPF saturate at this value irrespective of the 
number of particles used. 

B. Example 2 

A simple tracking system can be characterized by the 
O.OII3 and the sampling period 



following parameters, A = 
0.81 and P(0) 



T = 0.1. 

In Example 2, note that KLPF works with much fewer 
particles than in Example 1. One can attribute this to the 
much higher value of the optimal mean squared error in 
Example 1 than in Example 2, as can be seen from the plots. 

VI. Conclusions 

We propose a Kalman like particle filter (KLPF) that, 
in the examples studied, required moderately small number 
of particles and therefore can obtain close to optimal per- 
formance with a computational complexity comaparable to 
the conventional Kalman filter. An important open issue is 
to determine the number of particles necessary to closely 
approximate the optimal filter. As earlier observed in [12], 
the error covariance matrix of the optimal filter seems to 
follow the modified Riccati recursion introduced in [II]. 
Determining whether this is the case, and why, remains an 
interesting open question. 
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(a) SOI-KF, Alg 1 and KLPF 




...KLPF,r^^,, = 3 
— 5 level - KF 



(b) MLQ-KF, Alg 1 and KLPF 

Fig. 3. Both in (a) and (b), all the filters are close to optimal. 
KLPF achieves optimal performance with remarkably few particles 
and hence has a complexity of the same order as that of SOI-KF 
and MLQ-KF 
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Appendix 

A. The Generalized Closed Skew-Normal Distribution 

Definition 1: For x £ R", we define the generalized closed 
skew-normal distribution, GCSNk,n{x; fi, E, D, si, S2, A), as fol- 
lows 

GCSNk,n{x;fi,^,D,si,S2,A) ^ iV,- (a;; a*, E) 

"l>n {S1,S2;D{X - fi) ,A) 



Lfc,n(.) 



$„(si,s2;0,A + DEDT) 



(13) 



Nk{x; HjTi) is a fc-dim gaussian r.v and ^n(yi 



The sizes of the matrices involved follow accordingly. This is a very 
simple generalization of the closed skew-normal (CSN) distribution 
defined in [19]. Infact, it reduces to the CSN if yi — — oo 
in all its components. Naturally, it inherits most of its closure 
properties. Skew elliptical distributions generated a lot of interest 
because they provide a much needed tool to handle skewness 
in statistical modeling and have a good number of properties in 
common with the standard normal distribution, such as closure 
under marginlization and conditioning. We will briefly discuss some 
of the useful properties. 

Lemma 1.1: If X ~ Nk{n,T.), Z ~ iV„ (0, A) be 
independent, then X/ {D {X — jj.) — Z £ (si, S2)) ~ 
GC5iVfc,„(/i,E,D,si,S2,A) 

In the following proposition, we provide a stochastic character- 
ization of the GCSN which can be identified as a special case of 
theorem [3. II for state space systems. 

Proposition 1.2: Suppose 



< 



7/2"'' ) where 



V 



U 



■^Nk (0, E - ED^(A + DT.D'^) 
-Nn. (^si,s2,0,A + DED^ 



^DE| and 



are independent, where Np (a, b, 0, R) denotes a p-variate normal 
distribution truncated component- wise to lie in (a, b), where a and 
b are p-dim vectors specifying the boundaries of a cube in R'', 
i.e., Z^Np{a,b,0,R) =^ fz{z) = -^M|0^i^^^^_^j. 

Then, F = V + ED^ (A + DED^) f/ is distributed as 
GGSN (0,E,_D,i/, A). 

The proposition will be proved using lemmas 11.31 and 11.41 which 
are stated below. 

Lemma 1.3: If X ~ GG SNk,n(p,^, D , si, $2, A), then the 
moment generating function of X is given by 



Mxit) 



$„(si,S2;DEf, A + DED' 



1 T 

■ expl/it + —t Et) 



$„(si,s2;0,A + DEDT) 
Now, we will compute the moment generating function of a 
truncated multivariate normal distribution. 



Lemma 1.4: Let X ~ Np{a, b, 0, R), then the moment generat- 
ing function of X is given by 



Mxit) 



<in{si,S2;DT,t,R) 



exp 



r Rt 



$„(si,S2;0,i?) V 2 

We will now complete the proof of the proposition. 
Proof: [Proof of Proposition 11.21 



(14) 



Let S 



T,D' (A + DT.D' 



then Y 



V + SU and 



Mrit) = Mv{t)Msu{t). Now, Msu{t) = Mu{S'^t). Plugging 
in the formula for the moment generating function derived in lemma 
1.41 we see that My it) is the moment generating function of a 
GCS'iVfc,„(/i, S,L),si,S2, A) using lemma [73] ■ 
Proposition 1.5: Let X ~ GCSNk,n{fJ., S, L», si, S2, A), ~ 
Nk{l-i.w, Q) be independent and define Y = AX + W , then Y ~ 
GCSNk,nif^Y, Sy, Dy,si,S2, Ay), where 



Dy = s; 



and 



DYiD^ — Dy^y Dy , PY = M + M'v 



Ey = AT,A' + I 

Ay- = A 

This proof of the proposition follows directly from the following 
two lemmas. The first states that the GCSN is closed under full 
rank linear transformations while the second states its closure under 
addition of gaussian noise. 

Lemma 1.6: If Y ~ GCSNk,n{fJ.,T., D, si, S2, A) and A 
is a full column rank k x p {k > p) matrix, then 
AY ~ GCSNk,„{ApL,AT.A'^,D{A'^AyA'^,si,S2,A). If A 
is a full row rank p x k (p > k) matrix, then AY ~ 
GCSNk,„{fiA,T,A,DA,v, Aa), where 

fiA = A^i, Ea = AT.A^, Da = DEA'^T,^^ 

Aa = A + DED^ - DaT,aDa 
The following lemma states that the GCSN is closed under 
addition of gaussian noise. 

Lemma 1.7: If X ~ GCSNk,ni^l,T,, D, si, S2, A) and 
W ~ Nk{fiw,'Sw) are independent, then X + W ~ 
GCSNk,n {li-x+w, Ex+w, Dx+w, si, S2, Ax+w), 
where ij.x+w = A* + A'w, Ex+iv = E + Ew, 
-Dx+w = DE(E + Ew)"^ and Ax+w = A + DED"^ - 



i3 



'X+W'^X+W^'x + VK 



B. The Conditional State Density 

We will now show that the probability density of the state 
conditioned on the past quantized innovations / (a;(7i)/q„) follows 
a GCSN. Time update and Measurement update recursions will 
be provided for the parameters of the conditional density. All the 
parameters have very simple interpretations in terms of the statistics 
of the state space model. For ease of exposition, we will provide 
all the recursions for a time-invariant state space model. But the 
derivations can be trivially extended to the time varying case. 

Theorem 1.8 (Time and Measurement Updates): The 
probability density function of x{n)/q,n is given by 
GCSNk,m+i{x{n);0, E(?i/m), D{n/m), Si_„/„, S2,„/,„, 
A(n/m)), for m = n— 1, n and V n > 0. The recursions relating 
the parameters of the distributions of x{n— l)/Qn-i, x{n) / Q„-i 
and x{n)/Qn are given by 

E(n) = T,{n/n-l) = AT,{n~l/n~l)A'^ + W, 

D{n) = D{n/n - 1) = D(n-l/n-l)E(n-l/n-l)A^E(n)"^ 



given by the following recursions 

E(n/n) = E(n), D(n/n) 

, A{n/n) 



S2,n 



S2,n-1 
S2(n) 



Din/n-l 
H 

"A(n) 

af, 



Sl,. 



Sl,n-1 

Si{n) 



(16) 



where x(n) = x{n/n~l) = £'(a::(7i)/q„_i). In the rest of 
the appendix, we will use the notation R{n) — A(n/n) + 



D{n/n)T,{n)D{n/n) 
from Eq iflSll). 



T W 



A(n)-f D(n)E(n)D(n)^ ((*) follows 



C. Interpreting the Parameters 

Theorem 11.81 outlines how the conditional state density evolves 
with time. But, there is little intuition in the recursions proposed 
there. In the following lines, we will explain the connection of each 
of the parameters to the state space model, so that one can afford 
a better insight into the overall filtering technique first proposed in 
[9]. We will start by noting that E(n) is equal to the state covariance 
at time n. 

Lemma 1.9: E(n) = Ij.'c(n) 
The next Lemma states that R{n) — ||3^„|p, i.e., R{n) is the 
covariance of the measurement vector [j/(0), . . . , j/(n)]"^ 

Lemma 1.10: R(n) = ||3^„||^ = Ry{n) 
The following Lemma gives an expression for P(qn) 

Lemma 1.11: The probability of quantized innovations is given 

by 

P(q„) = $n+i(si, „, S2,„; 0, i?a(n)) 
In the following Lemma, we look at the cross covariance of x{n) 
and y„ 

Lemma 1.12: {yn,x{n)) = D{n/n)'^T}{n) 
Now, applying Proposition 1 1.21 on the conditional state density, we 
can write a'(n)/Q„=q„ stochastically as 

x(n)/q„ = Z + E(n)D(n/n)^(i?(n))"^yq(n), where 
Z - iVfe(0, E(?i) - T,{n)D{n/n){R{n))-'^ D{n/nfT,{n)) 
y^in) - iV„+i(si 0,Ry{n)) (17) 

Using Lemmas O [TTTol and WTtI one can see that 

E(n) - E{n)D{n/n)(R{n)y^ D{n/nfT}{n) = 

\\x{n)f - {x{n),y„) iWynfy' {yn,x{n)) and (18) 

j:{n)D(n/nf{R(n))-^ = {x(n),y„} [Wynf)'' (19) 

It is now easy to see that Eqs J18t and J19l > are consistent with 
theorem [3. II and hence we have 



xin)/q„ = Nk{0,Pkaiin/n))+R^y{n) {R{n))-^ y^{n) (20) 

where Pkai{n/n) is the error covariance of the Kalman filter at time 
step n. The methodology outlined in the appendix can be seen as 
an alternate way of arriving at Alg 3. However it is circumlocutions 
and not very intuitive. Moreover, the stochastic characterization in 
Eq l l20t cannot be extended to the more general case of theorem 
13.11 using the theory of GCSN. However, it serves as conclusive 
evidence that the conditional state density is not Gaussian, contrary 
to what is assumed in [9], [10]. 



_A 

S2,n-1 — S2,„/„_l — S2,„_l/„_l, 

A(n) = A(n/n-l) = A(n-l/n-l) + 

D(n-l/n-l)E(n-l/n-l)D(n-l/7i-l) - D^'E^.D^ (15) 
Eq l ll5t constitutes the time update. The measurement update is 



